R包EBSeq作差异基因分析

R包EBSeq差异基因分析

本篇简介R语言EBSeq包的差异基因分析流程。

EBSeq应用经验贝叶斯分层模型（empirical Bayes hierarchical model），寻找组间差异表达的gene（基因水平）或isoforms（可理解为由同一个基因编码的不同蛋白，由可变剪切产生）。

EBSeq简介：https://www.biostat.wisc.edu/~kendzior/EBSEQ/

安装EBSeq包

EBSeq可直接使用Bioconductor安装。

#Bioconductor 安装 EBSeq
BiocManager::install('EBSeq')

#加载
library(EBSeq)

后面有可能在加载EBSeq时会遇到这么一个问题，blockmodeling包载入时提示如下错误。

根据提示，我们在相应的路径下找到“CITATION”这个文件，然后用记事本等工具打开它，即可知道报错的原因：里面出现了R环境无法识别的字符（作者名称）。解决起来也很方便，把这些字符删掉就行了（删掉、或者替换成其它可以被R识别的字符都可以，而使用“#”注释掉是没用的）。不过由于无法被R环境识别的字符不止出现了一行，因此我们就简单粗暴点吧，把三个“citEntry()”里面的字符全删光得了（虽然觉得有些对不起作者……），肯定就没问题了。

EBSeq差异基因分析

下文以EBSeq包中自带的示例数据为例，简要展示差异基因分析流程。

值得一提的是，除了两组间的差异分析，EBSeq包中还提供了针对于三组及以上数据间的差异分析方法。以下分别以两组间差异基因分析以及三组间差异基因分析来演示（两组模式和多组模式的方法有些区别，详见下文）。

两组间差异基因分析

内置测试数据集GeneMat，是一个基因表达数据矩阵，包含10列样本共计1000行基因。我们载入该数据集，并手动将该数据集中的样本划分为两组，前5列的样本定义为C1组，后5列的样本定义为C2组，手动设置分组因子的顺序C1在前C2在后，在后续将比较C1中的基因相对于C2是上调了还是下调了。

#测试数据
data(GeneMat)

#手动将测试基因数据集设置为两组，前 5 列为 C1 组，后 5 列 C2 组
group <- factor(rep(c('C1', 'C2'), each = 5), levels = c('C1', 'C2'))

第一步，标准化因子

在差异基因分析前，需要分别对每个样本指定一个标准化因子（sample-specific normalization factor，在EBSeq中称为library size factor，简写为l_s），用于对数据标准化。标准化方法根据实际情况来确定，常用的标准化方法例如中位数标准化（Median Normalization）、分位数标准化（Quantile Normalization）、缩放标准化（Scaling Normalization）等。

此处对于示例数据基因表达数据矩阵GeneMat，选择使用了中位数标准化方法。结果返回一组向量。此处共包含10个元素（中位数标准化因子），分别对应了10个样本。

#通过 MedianNorm() 来获取标准化因子，等同于 DESeq 中的中位数标准化（Median Normalization）方法
ls_factor <- MedianNorm(GeneMat)

#其它标准化方法，例如通过 QuantileNorm() 实现上四分位数标准化
#ls_factor <- QuantileNorm(GeneMat, 0.75)
#等等方法，根据实际情况选择

第二步，差异基因分析

默认方法

之后使用EBTest()函数执行两组间的EBSeq差异分析。

#EBSeq 差异基因分析，本示例迭代 5 次是合适的，详见 ?EBTest
EBOut <- EBTest(Data = GeneMat, Conditions = group, sizeFactors = ls_factor, maxround = 5)
summary(EBOut)

#如果这些结果中，最后两个值相差小于 0.01 就表明迭代次数 maxround 合理
EBOut$Alpha
EBOut$Beta
EBOut$P

主要参数项：Data，指定基因表达矩阵GeneMat；Conditions，指定样本分组信息group，确保group中的分组名称顺序与GeneMat中的样本名称顺序（列顺序）对应，group为因子类型，因子顺序决定了比较的分组（本示例因子顺序，C1在前C2在后，在后续将比较C1中的基因相对于C2是上调了还是下调了）；sizeFactors，输入上述所得的样本标准化因子，以实现对数据的标准化；maxround = 5，对于本示例迭代5次计算，真实的数据分析中还需根据实际情况选择迭代次数。

3个指标用于判断选择的迭代次数maxround是否合理，一般来说最后两个值相差小于0.01就可以了。

FDR校正

通过GetDEResults()函数，获取差异基因。这里指定参数FDR = 0.05，意为为在执行FDR校正后以p<0.05为标准鉴定差异基因。

#鉴定差异基因（FDR < 0.05，其他参数如 Fold Change 等使用默认值），详见 ?GetDEResults
EBDERes <- GetDEResults(EBOut, FDR = 0.05)

结果赋值给列表EBDERes，包含3项内容。DEfound为鉴定的差异基因，对于示例数据我们共鉴定到了95个差异基因；PPMat包含两列信息，PPEE（within-condition variability的后验概率）和PPDE（across-condition variability的后验概率），如果PPEE<0.05（即可以理解为PPDE显著高于PPEE），即可将该基因定义为差异基因（DE gene）；Status则记录了每种基因的判别状态，包含DE（差异）、EE（无差异）、Filtered: Low Expression、Filtered: Fold Change、Filtered: Fold Change Ratio（后面三种表示该基因被过滤掉）五种类型。

既然根据PPEE<0.05的标准判别差异基因，那么即可知这里的PPEE其实就相当于FDR p-value了。也就是说上述方法中，我们默认根据FDR p-value < 0.05判断是否为差异基因。

#提取每个基因的 FDR p-value
FDR <- EBDERes$PPMat[ ,1]

提取Fold Change (FC)值

对于Fold Change (FC)值，可以通过PostFC()函数获得。

#获取获取每个基因的 Fold Change (FC) ，详见 ?PostFC
gene_FC <- PostFC(EBOut)

结果中包含RealFC和PostFC。二者的区别在于，RealFC基于原始数据获得，PostFC基于标准化后的数据获得，因此PostFC（posterior fold）会给低表达的基因提供一个较低的极值，以降低低表达基因的重要性。例如，如果gene1的mean1=5000且mean2=1000，则其RealFC和PostFC相差不大，大致均为5；如果gene2的mean1=5且mean2=1，则其RealFC为5，但PostFC将远低于5而接近1；此时当我们按PostFC进行排序时，gene2将不如gene1重要。根据实际需要，选择使用RealFC或PostFC。

可以通过PlotPostVsRawFC()比较每个基因的RealFC和PostFC。

#Plot PostFC vs RealFC
PlotPostVsRawFC(EBOut, gene_FC)

#不妨直接观察基因表达水平与 |RealFC - PostFC| 值的关系
FC_diff <- data.frame(C1_Mean = unlist(EBOut$C1Mean), C2_Mean = unlist(EBOut$C2Mean), FC_diff = abs(gene_FC$RealFC-gene_FC$PostFC))

下面的点图展示了示例数据中1000个基因的RealFC（FC）和PostFC（Posterior FC）的关系，基因的根据cross-condition mean（adjusted by the normalization factors）排序（升序排列，因此低表达的基因排在前面，展示为绿色；高表达基因排在后面，展示为红色）。越低表达的基因，RealFC和PostFC差别越明显，因此越偏离斜对角线。这就很明了地为我们展示了PostFC在缩小低表达基因差异方面的重要性。

或者，我们直接统计RealFC和PostFC的差值的绝对值，并结合各基因在C1组和C2组中的表达量一起展示。在结果中可以轻易地看出，越低表达的基因，RealFC和PostFC的差值的绝对值越大。

模型诊断

EBSeq依赖于参数假设，因此在每次分析时需要执行检查。可使用Q-Q图，拟合经验q's和使用Beta先验分布模拟的q's，评估参数假设Beta先验的合理性，如果两个条件下的数据点位于y = x线上，则表明Beta先验是合理的。或者，根据密度图来判断，如果二者趋势一致，则表明Beta先验合理。

#使用 Q-Q 图，拟合经验 q's 和使用 β 先验分布模拟的 q's，评估 β 先验合理性
par(mfrow = c(1, 2))
QQP(EBOut)

#查看经验 q's 和使用 β 先验分布模拟的 q's 的密度图，评估 β 先验合理性
par(mfrow = c(1, 2))
DenNHist(EBOut)

整合输出结果

确认以上结果没什么问题的话，不妨对以上的主要统计结果做个整合，输出在本地，便于后续自定义筛选等。

#整合版输出示例
C1_Mean <- unlist(EBOut$C1Mean) #C1 分组中，各基因标准化后的平均表达量
C2_Mean <- unlist(EBOut$C2Mean) #C2 分组中，各基因标准化后的平均表达量
FoldChange <- gene_FC$PostFC #使用 PostFC 获取
log2FoldChange <- log(FoldChange, 2) #log2 Fold Change

EBSeq_stat <- data.frame(C1_Mean, C2_Mean, FoldChange, log2FoldChange, FDR)
write.table(EBSeq_stat, 'EBSeq_stat_two.txt', sep = '\t', col.names = NA, quote = FALSE)

第一列为各基因的id；C1_Mean和C2_Mean列为两个分组中，各基因标准化后的平均表达值；FoldChange列就是Fold Change (FC)值，我们使用基于标准化后的数据获得的FC值，即PostFC；log2FoldChange就是对FoldChange作了log2转化；FDR，也就是上文中提到的PPEE值，这里相当于FDR p-value。

多组间差异基因分析

以三组间比较为例 

该R包同样提供了针对多组数据间的比较方法。这样，我们就用不着再单独执行两两分组间的比较并最后再整合了。

和上面的示例一样，这里同样使用EBSeq自带的数据集进行演示。测试数据集MultiGeneMat，是一个基因表达数据矩阵，包含6列样本共计500行基因。我们载入该数据集，并手动将该数据集中的样本划分为三组，前2列的样本定义为C1组，中间2列为C2组，最后2列为C3组。

#测试数据
data(MultiGeneMat)

#分组，模拟 3 组
group <- factor(rep(c('C1', 'C2', 'C3'), each = 2), levels = c('C1', 'C2', 'C3'))

第一步，标准化因子、预指定差异分析模式（类型）

首先同样是获取标准化因子（标准化因子的描述见上文）。

#通过 MedianNorm() 来获取标准化因子，等同于 DESeq 中的中位数标准化（Median Normalization）方法
ls_factor <- MedianNorm(MultiGeneMat)

多组间比较的方法与两组间比较有些区别，需要预先指定基差异分析模式。后面在计算时，对计算结果作判断，根据预先指定的差异分析模式划分归类。比方说本示例中共设置了3个分组（C1、C2、C3），那么基因差异模式就有5种。

#基因差异模式 Pattern
parti <- GetPatterns(group)
PlotPattern(parti)

Pattern1，某基因在3组间无差异；Pattern2，某基因在C3组与其它两组显著不同（无论是它高了还是低了，都算）；Pattern3，某基因在C2组与其它两组显著不同（无论是它高了还是低了，都算）；Pattern4，某基因在C1组与其它两组显著不同（无论是它高了还是低了，都算）；Pattern5，某基因在C1、C2、C3组之间均存在显著不同（无论该基因在各分组中的表达值C1>C2>C3、C1>C3>C2、C3>C1>C2等等，都算）。这个还是不难理解的吧。

在实际的分析中，我们可能不关注这么多情况，选择部分差异模式就可以了，并且全部用于计算也会相当消耗资源，那么也不妨过滤一些不必要的组合。比方说这里我们只想关注Pattern1和Pattern5这两类的话，可以这样筛选下。

第二步，差异基因分析

因为分组不多，数据也小，这里我们仍然分析了所有可能的差异类型组合。

默认方法

在大于两组的情况下，使用EBMultiTest()函数运行EBSeq差异分析流程。它的用法和上文中的EBTest()（两组情形）大致一致。

#EBSeq 计算，本示例迭代 5 次是合适的，详见 ?EBMultiTest
MultiOut <- EBMultiTest(Data = MultiGeneMat, Conditions = group, sizeFactors = ls_factor, AllParti = parti, maxround = 5)
summary(MultiOut)

主要参数：Data，指定基因表达数据矩阵MultiGeneMat；Conditions，指定样本分组信息group，确保group中的分组名称顺序与MultiGeneMat中的样本名称顺序（列顺序）对应，group为因子类型，但在多组比较时分组因子的顺序不重要；sizeFactors，输入上述所得的样本标准化因子，以实现对数据的标准化；AllParti，指定基因差异模式；maxround = 5，对于本示例迭代5次计算，真实的数据分析中还需根据实际情况选择迭代次数（迭代次数的选择是否合理，在上文的方法中已有提到）。

获取后验概率

在多组比较模式下，使用GetMultiPP()函数来获得每个基因在每种可能的差异模式中的后验概率。

#获得每个基因在每个 Pattern 中的后验概率，详见 ?GetMultiPP
MultiPP <- GetMultiPP(MultiOut)

结果赋值给列表MultiPP，包含3项内容。Patterns为先前预指定的基因差异模式；PP为各基因在各模式中的后验概率；MAP为最终划分的各基因的Pattern，根据PP作判断，一般来讲哪个Pattern的后验概率高就划分为哪种模式了。

模型诊断

可参考上文方法。

整合输出结果

我们对以上的主要统计结果做个整合，输出在本地，便于后续自定义筛选等。在多组比较模式下，各分组中各基因标准化后的平均表达量，会存储在同一个list中。对于上述示例，就是list“MultiOut$SPMean”，其中又包含了3个list [1]、[2]、[3]，分别对应了一开始指定的分组因子顺序（这里的顺序是C1、C2、C3），因此“MultiOut$SPMean”中的list [1]、[2]、[3]即分别对应分组C1、C2、C3。实际的数据分析中注意这点，不要提取错数据。

#整合并输出，一个示例
C1_Mean <- unlist(MultiOut$SPMean[[1]][[1]]) #C1 分组中，各基因标准化后的平均表达量
C2_Mean <- unlist(MultiOut$SPMean[[1]][[2]]) #C2 分组中，各基因标准化后的平均表达量
C3_Mean <- unlist(MultiOut$SPMean[[1]][[3]]) #C3 分组中，各基因标准化后的平均表达量

MAP <- MultiPP$MAP #Pattern 类型
EBSeq_stat <- data.frame(C1_Mean, C2_Mean, C3_Mean, MAP)
for (i in 1:nrow(EBSeq_stat)) EBSeq_stat[i,'C1_C2_C3'] <- paste(parti[EBSeq_stat[i,'MAP'], ], collapse = '_')

write.table(EBSeq_stat, 'EBSeq_stat_three.txt', sep = '\t', col.names = NA, quote = FALSE)

第一列为各基因的id；C1_Mean、C2_Mean、C3_Mean列为3个分组中，各基因标准化后的平均表达值；MAP和C1_C2_C3为Pattern类型。

可视化示例

ggplot2差异火山图（两组差异）

两组间差异的可视化，通常就是火山图了。如下展示了ggplot2绘制火山图的示例。

EBSeq_stat <- read.delim('EBSeq_stat_two.txt', sep = '\t')

#例如这里自定义根据 |log2FC| >= 1 和 FDR < 0.05 标记差异类型
EBSeq_stat[which(EBSeq_stat$FDR < 0.05 & EBSeq_stat$log2FoldChange <= -1),'sig'] <- 'down'
EBSeq_stat[which(EBSeq_stat$FDR < 0.05 & EBSeq_stat$log2FoldChange >= 1),'sig'] <- 'rich'
EBSeq_stat[which(EBSeq_stat$FDR >= 0.05 | abs(EBSeq_stat$log2FoldChange) < 1),'sig'] <- 'no diff'

#横轴 log2FC（使用 PostFC 获取），纵轴 log10 转化后的标准化后基因表达量的平均值，颜色表示差异
library(ggplot2)

volcano_plot <- ggplot(EBSeq_stat, aes(log2FoldChange, log10((C1_Mean + C2_Mean)/2))) +
geom_point(aes(color = sig), alpha = 0.6, size = 1) +
scale_colour_manual(values = c('blue2', 'red2', 'gray'), limits = unique(EBSeq_stat$sig)) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.25) +
geom_hline(yintercept = -log(0.05, 10), color = 'gray', size = 0.5) +
xlim(-3, 3) +
labs(x = 'log2 Fold Change (C1 vs C2)', y = 'log10 Mean of Normalized Count', color = '')

#ggsave('volcano_plot.pdf', volcano_plot, width = 5, height = 5)
ggsave('volcano_plot.png', volcano_plot, width = 5, height = 5)

ggtern三元相图（三组差异）

三组间差异的可视化，可以考虑使用三元相图呈现。如下展示了ggtern（ggplot2的衍生包，语法和ggplot2几乎一致，如果你也喜欢ggplot2的风格，就可以考虑它）三元相图示例。

EBSeq_stat <- read.delim('EBSeq_stat_three.txt', sep = '\t')

#由于多组模式无 Fold Change (FC)，因此我们直接根据各基因在各组中的平均表达量判别组间表达水平高低
#同时结合 Pattern，判别显著性

#计算各基因在哪个分组中的表达最高
EBSeq_stat$max <- apply(EBSeq_stat[2:4], 1, function(x) names(x)[x %in% max(x)])

#根据 Pattern 以及各基因所在表达最高的分组，标记基因富集分组
EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern1'),'high_sample'] <- 'none'

EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern2' & EBSeq_stat$max == 'C3_Mean'),'high_sample'] <- 'C3'
EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern2' & EBSeq_stat$max %in% c('C1_Mean', 'C2_Mean')),'high_sample'] <- 'C1 & C2'

EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern3' & EBSeq_stat$max == 'C2_Mean'),'high_sample'] <- 'C2'
EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern3' & EBSeq_stat$max %in% c('C1_Mean', 'C3_Mean')),'high_sample'] <- 'C1 & C3'

EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern4' & EBSeq_stat$max == 'C1_Mean'),'high_sample'] <- 'C1'
EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern4' & EBSeq_stat$max %in% c('C2_Mean', 'C3_Mean')),'high_sample'] <- 'C2 & C3'

EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern5' & EBSeq_stat$max == 'C1_Mean'),'high_sample'] <- 'C1'
EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern5' & EBSeq_stat$max == 'C2_Mean'),'high_sample'] <- 'C2'
EBSeq_stat[which(EBSeq_stat$MAP == 'Pattern5' & EBSeq_stat$max == 'C3_Mean'),'high_sample'] <- 'C3'

#这里计算了各基因在所有样本中表达量的均值的平方根，用于作图时定义点的大小......
EBSeq_stat$size <- apply(EBSeq_stat[2:4], 1, function(x) sqrt(mean(x)))

#ggtern 三元相图示例
library(ggtern)

tern_plot <- ggtern(EBSeq_stat, aes(C1_Mean, C2_Mean, C3_Mean)) +
geom_mask() +
geom_point(aes(color = high_sample, size = size), alpha = 0.6) +
scale_size(range = c(0, 4)) +
scale_colour_manual(values = c('red', 'blue', 'green3', 'purple', 'orange', 'skyblue', 'gray'), limits = c('C1', 'C2', 'C3', 'C1 & C2', 'C1 & C3', 'C2 & C3', 'none')) +
theme_bw() +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
labs(x = 'C1', y = 'C2', z = 'C3', color = 'en-riched') +
guides(size = FALSE)

#ggsave('tern_plot.pdf', tern_plot, width = 6.5, height = 5)
ggsave('tern_plot.png', tern_plot, width = 6.5, height = 5)

友情链接